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ABSTRACT 


A  method  has  been  developed  to  investigate  the  sensitivity  of  the  solutions 
of  large  sets  of  coupled,  non-linear  rate  equations  to  uncertainties  in  the  rate 
coefficients.  This  method  is  based  on  varying  all  the  rate  coefficients  simul¬ 
taneously  through  the  introduction  of  a  parameter  in  such  a  way  that  the  output 
concentrations  become  periodic  functions  of  this  parameter  at  any  given  time  t. 
The  concentration  of  the  chemical  species  are  then  Fourier  analysed  at  time  t. 

We  show  via  an  application  of  Weyl's  ergonc  theorem  that  a  subset  of  the  Fourier 

/dcA 

coefficients  is  related  to  h  the  rate  of  change  of  the  concentration  of 
species  i  with  respect  to  the  rate  constant  for  reaction  I  averaged  over  the 
uncertainties  of  all  tne  other  rate  coefficients .  Thus  a  large  Fourier  coeffi¬ 
cient  corresponds  to  a  large  sensitivity,  a  small  Fourier  coefficient  corresponds 
to  a  small  sensitivity.  The  amount  of  numerical  integration  required  to  calcu¬ 
late  these  Fourier  coefficients  is  considerably  less  than  that  required  in  test: 
of  sensitivity  where  one  varies  one  rate  coefficient  at  a  time,  while  homing 
all  others  fixed.  The  Fourier  method  developed  in  this  paper  is  not  limited 
to  chemical  rate  equation,  but  can  be  applied  to  the  study  of  the  sensitivity 
of  any  large  system  of  coupled,  non-linear  differential  equations  with  respect 
to  the  uncertainties  in  the  modeling  parameters. 


I.  introduction 


Sets  of  coupled,  non-linear  rate  equations  arise  in  a  number  of  disciplines. 

A  classic  example  is  t ha t  of  chemical  rate  equations.  In  the  study  of  combustion, 
air-pollution,  upper  atmosphere  phenomena  and  chemical  lasers  as  many  as  100 
coupled  rate  equations  involving  some  50  separate  species  may  be  needed  to 
account  for  the  properties  of  such  systems.  One  is  then  faced  with  the  problem 
of  solving  a  large  set  of  coupled,  non-linear  differential  equations  of  the  form 


dc  * 
3t  =  F 


c .  { k  H 


(l.D 


when  c  is  a  vector  of  concentrations,  { k}  a  set  of  rate  coefficients  and  P  some 
complicated  function  of  the  concentrations.  While  one  may  prove  existence  and 
stability  theorems  for  the  equilibrium  point^,  tne  only  way  to  solve  these 
equations,  i.e.  to  obtain  all  the  species  concentrations  as  functions  of  time, 
is  through  the  use  of  a  high-speed  computer. 

Unfortunately,  the  rate  coefficients  {or  cross  sections)  for  many  reactions 
are  not  known  with  high  accuracy  and  indeed  may  be  uncertain  by  one  or  more 
orders  of  magnitude.  This  gives  rise  to  the  very  important  problem  of 
"sensitivity"  wtiich  may  be  defined  as  the  effect  of  uncertainties  in  the  rate 
coefficients  on  tne  calculated  concentrations  of  all  the  various  intermediate 
and  product  species.  The  uncertainty  in  the  rate  coefficients  of  certain 
"important"  reactions  in  the  reaction  scheme  may  have  a  significant  effect  on 
tne  output  functior  (for  instance,  concentration  at  time  t),  while  uncertainties 
of  the  same  magnitude  in  rate  coefficients  of  "unimportant"  reaction  in  the 
reaction  seneme  may  hardly  effect  the  output  function.  The  reliability  o  the 
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output  numbers  clearly  cannot  be  established  without  knowledge  of  the  sensitivity 
of  the  output  data  to  the  uncertainty  in  the  input  parameters. 

The  problem  is  to  find  a  practical  method  of  determining  the  effects  of 
the  uncertainties  in  the  rate  coefficients  on  the  solutions  of  the  rate  equations. 
Since  we  are  interested  in  the  situation  where  the  uncertainties  in  tne  rate 
coefficient  may  be  orders  of  magnitude,  linearisation  schemes  are  not  appropriate. 
A  "brute  force"  method  of  investigating  the  sensitivity  is  not  feasible  as  can 
readily  be  seen  from  the  following  example.  Suppose  we  have  a  reaction  scheme 
which  has  n  coupled  reactions  involving  m  different  chemical  species.  Let  us 
furthermore  assume  that  we  wish  to  calculate  the  concentrations  of  the  m  species 
at  some  time  t  for  z  different  values  each  of  the  2n  rate  coefficients.  If  we 
now  change  one  rate  coefficient  at  a  time  while  holding  all  the  others  fixed, 
we  would  lave  to  carry  out  z^n  integrations  of  the  rate  equation  (1.1)  to  time  t. 
For  the  m  different  species,  this  procedure  will  give  rise  to  a  pnnt-out  of 
m(z)^n  concentrations.  If  we  know  to  a  good  accuracy  the  equilibrium  constants 
for  all  the  reactions  and  apply  the  principle  of  detailed  balance,  the  number 
of  independent  rate  coefficients  will  be  reduced  to  n  and  we  would  have  to  carry 
out  zn  integrations  up  to  time  t  for  eacn  species  m.  In  either  case,  for  n 
large,  it  is  obvious  that  the  time  and  expense  involved  in  such  an  analysis  of 
sensitivity  is  prohibitive  and  the  print-outs  so  numerous  as  to  defy  a  simple 
analysis  of  the  results.  Clearly,  one  needs  to  devise  some  more  powerful 
method  for  the  study  of  sensitivity. 

Our  approach  to  this  problem  is  to  ask  for  a  less  detailed  description  of 
the  effect  of  rate  coefficient  uncertainty  on  the  output  function  at  any  given 
time.  He  vary  all  the  rate  coefficients  simultaneously  so  as  to  explore  the 
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entire  space  of  uncertainties  in  the  rate  coefficient  set  { k}  .  As  will  be  seen 
below,  this  turns  out  to  be  equivalent  to  varying  a  single  rate  coefficient  and 
then  averaging  the  attendant  concentration  changes  over  the  uncertainties  of  all 
the  other  rate  coefficients,  wnere  we  express  this  uncertainty  in  terms  of  a 
probability  distribution.  Our  approach  is  thus  related  to  a  "mean  field"  theory 
where  we  represent  the  fluctuations  of  the  field  by  the  uncertainties  in  the 
(n-1)  rate  coefficients  over  which  we  average. 

To  carry  out  this  program,  we  relate  eacn  rate  coefficient  k.  to  a  frequency 
um  and  introduce  a  parameter  s  which  simultaneously  varies  all  the  rate  coef¬ 
ficients  in  such  a  way  that  the  concentrations  a_t  a  fixed  time  become  periodic 
functions  of  s.  The  concentrations  can  then  be  Fourier  analyzed.  We  then  show 
that  a  certain  subset  of  these  Fourier  coefficients  can  be  related  to  the  first 
partial  derivative  of  tne  concentration  c.  of  species  i  with  respect  to  a  rate 
coefficient  k  averayed  over  the  uncertainties  of  all  tne  other  rate  coefficients. 

A  large  value  of  the  Fourier  coefficient  A'1'  tnen  shows  that  (~~-)is  large, 

“*£  \dk;/ 

i.e.  the  effect  of  i  change  in  the  £'th  rate  coefficient  on  the  concentration 
of  species  i  is  significant.  Conversely,  a  small  Fourier  coefficient  A^ 


indicates  tnat^-^~y  is  small,  i.e.  the  effect  of  the  variations  of  the  j'th 
rate  coefficient  on  the  concentration  of  species  i  is  small.  In  order  to  calcu¬ 
late  these  Fourier  coefficients,  we  must  integrate  the  rate  equations  numerically 
up  to  the  desired  time  for  each  value  of  the  parameter  s.  The  number  of  s  values 
which  we  include  in  our  parameter  set  determines  the  accuracy  to  which  we  can 
calculate  the  Fourier  coefficient;  the  larger  the  set  of  s  values,  the  more 
accurate  the  determination  of  the  Fourier  coefficients. 

Since  we  must  still  perform  numerical  integrations  of  the  rate  equations 


dCi\  ... 
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(1.1),  the  question  arises  wny  this  method  is  to  he  preferred  to  the  more 
direct  method  of  varying  each  rate  coefficient  separately  while  Keeping  all 
others  fixed.  As  will  be  shown  in  paper  II,  which  deals  primarily  with  the 
computer  calculations,  trie  number  of  integrations  required  by  the  Fourier 
method  is  0(nr)  where  r  is  a  small  integer  (r<10)  which  depends  upon  the  choice 
of  tne  frequencies  i  =  l,2,...,n.  It  can  readily  be  verified  tnat  for  n, 
z  large,  nr<(z)n.  The  computational  economy  of  this  method  of  analysis  thus 
becomes  more  pronounced  as  n,  the  number  of  reactions,  increases.  The  reason 

for  this  reduction  in  the  number  of  required  integrations  up  to  time  t  is  to 
be  found  in  the  fact  that  in  the  Fourier  method  we  sample  the  c(k)  space  at 
a  set  if  points  determined  by  the  values  of  the  set  { s }  and  the  vector  a), 
whereas  the  "brute  force"  method  involves  neither  sampling  nor  the  simultaneous 
variation  of  the  set  of  rate  coefficients  ik)  and  thus  requires  many 

more  integrations  of  the  rate  equations.  As  will  be  clear  from  the  body  of 
the  paper,  this  sampling  in  a  certain  sense  corresponds  to  the  averaging  over 
all  the  rate  coefficients.  The  reduction  in  the  number  of  required  integrations 
is  thus  intin, ately  related  to  the  use  of  a  "mean  field"  theory. 

Our  results  are  presented  in  two  papers,  I  and  II.  In  paper  I,  we  present 
the  theoretical  basis  of  our  method  without  explicit  reference  to  the  verifying 
computer  experiments.  In  paper  II,  we  present  the  results  of  our  computer 
calculations  as  well  as  a  detailed  discussion  of  the  problems  involved  in  such 
cal culations . 

It  should  be  pointed  out  that  the  utility  of  this  Fourier  analysis  method 
of  testing  sensitivity  extends  beyond  the  confines  of  chemical  kinetics  and 
beyond  the  confines  of  differential  equations.  Large  sets  of  coupled,  non¬ 
linear  equations  are  used  in  many  fields  such  as  economics,  population  dynamics, 
weather  forecasting,  systems  analysis,  operations  research,  etc.  for  modeling 


and  predictive  purposes.  As  has  been  pointed  out  by  a  number  of  investigate 
it  is  important  that  sensitivity  tests  be  carried  out  on  such  systems  to 
identify  the  critical  parameters  and  to  validate  the  applicability  of  the 
models.  The  method  developed  iiere  can  be  applied  to  any  set  of  equations, 
differential,  integral,  algebraic,  etc.,  which  yield  an  output  as  a  complex 
finction  of  many  parameters.  In  fact,  we  have  used  an  analytic  function  of 
several  variables  to  test  some  of  our  ideas.  We  plan  to  apply  this  Fourier 
method  to  the  analysis  of  other  complex  systems  in  the  near  future. 


II.  Fourier  Analysis 


The  rate  equation  for  coupled  chemical  reactions  which  obey  the  law  of 
mass  action  can  be  written  in  the  form 


dc . 

i 

ar  = 


r,1 


(2.1) 


H  I 

where  c -(k, .,k  :t)  is  the  concentration  of  species  i  at  time  t,  v  •  =  v  •  -  v  • 

is  the  stoichiometric  coefficient  of  the  species  i  in  reaction  r,  wi  th  r=l,...,n 

labeling  the  different  reactions  in  the  reaction  system, and  where  k  (k  )  is 

r  -r 

I 

the  forward  (backward)  rate  coefficient  for  reaction  r.  The  coefficients  v  . 

It 

and  v r  j  are  non-negative  integers  so  that  the  stoichiometric  coefficient 
vr-  defined  above  can  be  positive,  negative  or  zero.  From  the  form  of  the  rate 
equations  (2.1)  one  can  prove  that  the  concentrations  are  bounded  and  non-negative, 
that  they  are  continuous  functions  of  the  time  and  the  rate  coefficients  and 
that  an  equilibrium  point  exists.  ^ 

In  order  to  determine  the  effects  of  uncertainties  in  the  rate  coefficients 
k+r  on  the  concentrations  c-  at  time  t  a  systematic  method  for  varying  the  k's 
must  be  developed.  We  define 

^  =  kj°>eui  0  =  1,. ...n)  (2.2) 

and 


Uj  =  (sin  i^s) 

where  k|°^is  the  "best  value"  of  the  rate  coefficient  (i.e. 
investigator  oelieves  to  be  the  best  available  value  based 
calculations),  the  "frequency"  w .  is  a  positive  integer,  s 


(2.3) 

the  one  which  the 
on  experiments  or 
is  a  parameter  and  f. 


-7- 


is  a  function  to  be  determined.  The  introduction  of  the  parameter  u.  and  the 
form  of  Eq.  (2.2)  permits  one.  to  effect  readily  order  of  magnitude  changes  in 
the  rate  coefficients.  The  form  of  equation  (2.3)  permits  one  to  vary  simul¬ 
taneously  all  the  rate  coefficients  by  varying  the  parameter  s.  The  rapidity 
of  the  variation  is  determined  for  each  by  the  magnitude  of  the  w-'s.  The 
u>  •  are  chosen  to  be  positive  integers  in  order  that  the  concentrations  at  a 
fixed  time  t  become  periodic  functions  of  s  with  period  2n: 

c  n-  ( s ;  t )  =  ci  (s  +  2n  ;t.)  (2.4) 

In  the  development  to  follow  we  suppress  this  dependence  of  the  concentrations 
on  time,  it  beiig  understood  that  the  analysis  is  carried  out  for  a  fixed  time  t. 

The  concentration  as  a  function  of  s  describes  a  closed  path.  That  is, 
for  eacn  i  and  for  every  value  of  s  we  oDtair,  a  point  in  k  space  with  value 
c^(l<);  as  s  changes  by  2n  we  return  to  the  same  values  of  t:  and  u,  [see  eqs. 

(2.2)  and  (2.3)],  and,  from  equation  (2.4),  to  the  same  value  c .  ( s ) .  Since 
c^(s)  is  periodic  on  2tt  wo  may  expand  it  in  a  Fourier  series 

A(  i)  ” 

ci(s)  =  -r-  +  S  (a 

1  *  r=l  V 

In  the  analysis  presented  below  we  are  interested  only  in  the  Fourier  sine 
coefficients  that  correspond  to  the  original  input  frequencies  uk ,  i.e.  the 

coefficients  given  by 


(i) 

r 


sin  rs  +  b 


*0) 


cos  rs 


(2 


5) 


“t  n 


f 


^  ^  =  -  f  c.(s)  sin  u  s  ds 


(2.6) 


i  =  1 ,2, . . . ,n 


We  now  wish  to  relate  tne  Fourier  coefficients  A'  to  the  effect  of  the 

u) . 


(i) 
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variation  of  tne  rate  coefficients  k^  and  the  concentrations  c..  In  order  to 
do  tnis  it  is  necessary  to  relate  the  s-  space  integral  of  Eq.  (2.6)  to  a  n- 
d  ir.iens  ional  integral  over  the  entire  u  space.  It  can  be  shown  that  if  the 
entire  u  space  is  covered  densely  as  s  is  varied,  then  the  integrals  over  tne 
s-  space  and  the  u  space  yield  identical  results,  fnis  is  just  a  statenent 
of  tne  ergodic  tneorem  which  permits  one  to  equate  time  and  phase  space  averages 
in  statistical  mechanics.^  If  the  frequencies  in  Eq.  (2.3)  were  chosen  to 
be  incommensurate,  then  the  function  c.(s)  i.e.  the  concentration  as  a  function 
of  the  parameter  s,  would  be  an  almost  periodic  function.  This  implies  that 
the  path  in  the  space  Rn(4),  where  the  0 1 s  are  defined  by 

0 .  =  ij.s  (mod  2")  (2.7) 

returns  arbitrarily  close  to  any  initial  point  as  s- — *<*>.  One  could  then  equate 
the  s-  space  integral  with  the  integral  over  the  n-  dimensional  3  space  as  was 
first  proved  by  Weyl . ^  However,  the  use  of  an  incommensurate  set  of  frequencies 
u  in  Eq.  (2.3)  would  require  that  the  numerical  evaluation  of  the  Fourier 
coefficients  A  of  Eq.  (2.6)  be  carried  out  over  an  infinite  period.  This 
clearly  is  not  a  feasible  procedure  on  any  known  computer. 

It  is  for  this  reason  that  we  introduce  integei  frequencies  ^  in  Eq.  (2.3). 
Their  use  leads  to  a  finite  period  (0  to  2«)  analysis  which  can  readily  be 
handled  on  a  computer.  Tnis  will  give  rise  to  an  error  in  the  analysis  since 
now  the  "phase  point"  will  no  longer  densely  cover  the  u  space  as  s  is  'aried 
and  the  s-space  and  u  space  integrals  therefore  do  not  yield  identical  results. 


Let  us  for  tne  moment  ignore  this  error  and  use  Weyl's  theorem  for  our 


periodic  function  c(s).  We  write 


Ai)  _  2 


/.2n  /-2tt  n 

—  r  ...  n  -v, s,n 

Jo  J0  j  =  l 


(2.8) 
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This  expression  can  be  put  in  a  mere  suggestive  fori  as  follows.  The  i ' th 
integral  is  rewritten  by  integration  by  parts  as 


1 

2n 


d0t 


.0^}  sin  0 


£ 


1_ 

2n 


d0 


dc_ 

l  S0„ 


cos 


(2.9) 


where  the  boundary  term  has  vanished  by  the  periodicity  of  the  integrand.  The 
use  of  Lq.(2.9)then  permits  us  to  rewrite  Eq.(2.8)  as 

•2t  r2i  ..  be- 


k(1)  .  _1 


u) 


(2tt  ) 


■J' 

•'0  •'o  J =  1 


d0. 


J  bQ 


cos  0 


(2.10) 


Tne  6  space  integral  must  now  be  related  to  the  u  space  integral.  The  transfor¬ 
mation  is  (see  Eqs.  (2.3'  and  (2.7)) 


du. 


bf -(sin  0j) 


cos  0 .  d0 . 
*3  J 


(2.11) 


where  we  require  that  f.  be  a  mono  tonic  function  of  its  argument.  Then 

J 


fM+D  ,f„ 

=  ...f 

a)  n  J  I 

1  ’  f,  (-i )  f„(-n 


'f"H)  dUj  cos  ot 

bf^sin  0^) 


(2.12) 


n'  ' '  [\  cos  0^ 


j  b  sin  9j 


In  order  to  obtain  a  definite  relation  between  and  bc-/bu0  we  must  choose 

1  *- 

some  particular  form  for  tne  function  fj (sin  Oj ) .  There  are  several  choices  which 
lead  to  useful  expressions  for  the  A^;  as  will  be  seen  below,  a  particularly 
advantageous  choice  is  to  set 
afjtsin  ep 


b  sin  0 


J 


2n  1 
cos  0.  =  — 
1  a  ■ 
J  J 


(2.13) 


where  a  is  a  parameter.  The  use  of  Eq.  (2.13)  in  Eq .  (2.12)  leads  to 
J 
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,(i>  =  2 


>.4*  n  ac. 

7  J  •  j  TT  duj  aj cos  °j  ■ 

•'-on  *-oo  J  =  | 


(2.14) 


Integration  of  Eq.  (2.13)  yields 


i  fj  +  sin  0.1 

fj(s1n  V  ■  uj =  w:  l"  r^rnnrj 

J  J 


(2.15) 


and  also  indicates  the  range  of  integration  in  Eq.  (2.14).  Expressing  cos  9^ 
as  a  function  of  Uj  from  Eq.  (2.15)  gives  for  Eq.  (2.14) 


A«>=-U  ...I  IT  du 

w £.  a.7rn  ■'-«>  *^-on  i  =  1  ^ 


cosh  a.u.  au 
J  J  *- 


(2.16) 


Since 


/-4°  n  a- 

I  ‘  '  I  TT  duj  cosh  a’TTTT 

-'-CO  •'~X>  j=l  1  1 


we  obtain  as  our  final  result 


Oc./du  > 

UJl  a£  1  Z 

where  the  bracket  in  E.q.  (2.18)  is  defined  by 

J  •j  7\  P<uj;  aj)  V(u, - ,un)  du, 

<»<", . - — 

...  I  TT  p( u . ;  a .)  du. 

J. L™>  -L*.  -;=]  J  J  1 


(2.17) 


(2.18) 


(2.19) 


p(uj:'aj)  cosh  a^u. 


(2.20) 


The  function  p(u.;a.)  can  be  interpreted  as  a  distribution  function  in  u  space 
J  J 

which  weights  the  uncertainty  in  the  rate  coefficients.  Equation  (2.18)  is  the 


-11- 


desired  relation  between  the  Fourier  coefficient  A^  and  the  chanqe  in  the 

wJi, 

concentration  of  species  i  with  a  change  in  the  Jt'th  rate  coefficient,  dc./du  , 

1  A* 

averaged  over  the  changes  of  all  the  other  rate  coefficients. 

The  particular  form  of  f .  (sin  9.)  that  we  have  chosen,  Eq.  (2.15)  has 
lead  to  a  weight  function  p(u.;a.)  for  each  rate  coefficient,  Eq.  (2.20),  which 

J  J 

has  several  convenient  properties.  As  a  function  of  u.,  the  function  p(u.;a.) 

J  J  J 

is  symmetrically  bell  shaped  about  u .=0  corresponding  to  k.=k(0^,  the  "best" 

J  J  J 

value  of  the  rate  coefficient  k..  As  a.,  which  is  a  parameter  at  our  disposal, 

is  increased,  the  weight  function  p(u.;a.)  narrows  about  u.=0;  this  corresponds 

J  J  J 

to  a  decreased  spread  in  the  values  of  the  rate  constant  k.,  i.e.  it  corresponds 

J 

to  a  narrower  range  of  the  uncertainty  of  k.  around  k(°).  In  the  limit  as  a 

J  J  j 

one  obtains 


lim  p(u.;a.)  =  n6(u.) 

a-. — J  J  J 

J 


(2.21) 


where  6(x)  is  the  delta  function  of  argument  x  which  implies  that  the  rate  con¬ 
stant  is  known  with  certainty  to  be  k^0^.  When  some  information  is  available 

on  the  spread  and  distribution  about  the  value  of  the  rate  coefficient  k.  about 

/  \  J 

its  "best"  value  ky',  one  can  determine  the  a,  through  the  standard  deviation 

J  J 

of  the  values  in  u  space 


(2.22) 


(2.23) 
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The  parameter  thus  permits  us  to  introduce  explicitly  the  spread  of  uncertain¬ 
ties  in  the  values  of  the  rate  coefficients  into  our  analysis. 

To  calculate  the  Fourier  coefficients  we  must  first  choose  a  suitable 

(ij 

i 

set  of  Uj.'s  and  a.'s.  What  constitutes  a  “suitable"  set  ofu's  will  be  discussed 
in  the  next  section  in  conjunction  with  the  error  analysis  and  in  more  detail 
in  paper  II.  The  set  of  u^'s  defines  the  u/s  as  a  function  of  s  according  to 
Eqs.  (2.3),  (2.7)  and  (2.15).  For  each  value  of  s  one  obtains  a  value  of  the 
concentration  ci  (k] , . . .  ,kn) .  In  principle  one  can  then  compute  the  to  any 
desired  accuracy  by  taking  enough  values  of  s, 


.m. 


(2.24) 


where  m  is  some  integer. 

It  is  important  to  point  out  that  our  main  interest  is  in  identifying  those 
rate  coefficients  whose  variation  significantly  effects  the  concentration  c. 
of  a  species  i  at  time  t,  and  those  rate  coefficients  whose  variation  has  only 
a  minimal  effect  on  the  species  concentration  c^ .  Thus  if  one  of  the  Fourier 
coefficients,  say  A^,  is  one  or  more  orders  of  magnitude  larger  than  all 

other  coefficients  A^,  i  =  1 ,2, . . .  ,n,  (^j),  than  the  variation  of  the  j'th 

i 

rate  coefficient  k j ,  clearly  has  a  larger  effect  on  the  concentration  c.(t) 
than  the  variation  of  the  other  rate  coefficients.  If  on  the  other  hand  all 

the  coefficients  A^  ^ ,  i  =  l,2,...,n,  are  of  the  same  order  or  magnitude, 

i 

then  the  concentration  of  species  i  at  time  t,  c.(t),  is  effected  essentially 
equally  by  the  variation  of  any  of  the  rate  coefficients  k  . 

One  problem  with  the  above  analysis  must  be  pointed  out.  The  Fourier 

coefficient  A  ^  of  Eq.  (2.18)  may  be  small  either  because  be. /du  is  small  over 

a  1  1 
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. ...  ^  i * 


the  whole  range  of  integration  of  Eq.  (2.19),  which  is  the  case  discussed  atuvj, 
or  because  dc^/du  chai  ges  sign  one  or  more  times  in  the  range  of  integration. 

Thus,  a  small  value  of  A  ^  does  not  necessarily  imply  that  the  concentration  c. 
is  insensitive  to  changes  in  the  rate  coefficients.  The  remedy  for  this  ditn- 
culty  would  be  to  find  a  relation  between  the  Fourier  coefficients,  or  a.  combination 
of  Fourier  coefficients,  and  an  everywhere  positive  (or  everywhere  negative) 
function  of  the  rate  of  change  of  the  concentrations  with  rate  coefficients,  such 
as  for  instance  ((dc-/du  )*").  We  have,  however,  not  been  able  to  establish 
such  a  relationship  and  it  seems  doubtful  that  a  simple  relationship  of  this  form 
exists.  In  carrying  out  computer  calculations  it  should,  however,  not  be  too 
difficult  to  verify  whether  dc^du^  at  any  given  time  t  is  monotone  or  not  in 

the  range  of  integration  over  the  u  space. 

Thus  while  one  can  assert  that  a  large  implies  high  sensitivity  of  the 

concentration  of  species  i  with  respect  to  changes  in  the  rate  coefficient  k^, 
the  converse  statement  does  not  necessarily  fol loty  wi thout  checking  on  the  mono¬ 
tonicity  of  dc^/du  as  discussed  above. 
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III.  Choice  of  the  Frequencies  w  •  • 


As  we  have  Mentioned  aoove,  the  application  of  Weyl's  theorem  in  going  from  s 
space,  Eq.  (2. fa),  to  the  u  space,  Eq.  (2.16),  via  the  6  space,  Eq.  (2.8),  must  iPad  to 
an  error  in  the  analysis  since  we  use  commensurate  (integer)  frequencies.  Tnis 
error  can  be  minimized  by  a  judicious  cnoice  of  the  integer  frequencies  and 
the  number  and  magnitude  ot  the  n  values  of  the  parameter  s  given  by  tq.  (2.24). 

We  limit  ourselves  here  to  some  qualitative  remarks  which  will,  however,  clearly 

indicate  the  nature  of  tne  pro  Mem. 

The  integer  frequencies  lead,  according  to  Eqs.  (2.7)  and  (2.24,, 
covering  of  the  n  dimensional  6  space  by  an  array  of  points  as  q  takes  on  its 
integer  values  ,=1.2... J».  Clearly  one  obtains  a  better  coverage  of  tne  t>  space 
and  thus  reduces  the  error  in  applying  Weyl's  theorem  if  one  can  increase  the 
density  of  points  in  the  n-dimensional  nypercute  and  if  one  can  distribute  the 
points  uniformly  within  the  space.  Since  2  and  q  are  both  integer,  the  points 
J. s  (mod  2,)  will  form  a  regular  lattice  in  i  space.  Our  objective  will  be  to 
make  this  point  lattice,  which  is  completely  generated  by  a  unit  cel  1,  as  uniform 
as  possible  in  all  n  diiections  by  a  judicious  choice  of  the  »  s.  As  our  measure 
of  uniformity  we  take  a.  nyfiercuoicjnU^li.  Without  further  information  about 
tne  behavior  of  the  output  function  this  seems  the  most  reasonable  choice. 

For  a  fixed  number  m  of  s  points,  chosen  for  computational  convenience, 
we  will  try  to  find  a  vector  2  (of  tne  infinite  set  of  2's)  which  gives  rise 
to  a  nypercubic  unit  cell.  Once  having  done  this  we  can  assert,  without  loss 


of  generality,  that  &  lies  along  one  edge  of  tins  hypercubic  unit  cell  whose 


length 


in 

2^[2L  There  are  now  two  ways  to  compute  tne  volume  V  of  the 

1  ^  m 
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■kb. 


cell.  In  the  n  dimensional  6  space,  it  follows  front  the  above  cons tr't. non 


\  m  f 


(3.1) 


3-t  wt  also  know  that  the  total  volume  of  the  n-dimensiona!  5  space  is  (2  )n. 
Since  there  are  m  unit  cells  in  that  space,  it  follows  that 


V  -  ( 2  n ) ri 
rn 


(3.2) 


Equating  (3  1)  and  (3.2)  then  yields 


I  u  I  -  l.J 


(3.3) 


tor  t n 2  relation  between  the  minimum  length  of  tne  vecto'  j  ,  the  dimension  n 
and  tne  number  of  s  points,  m,  for  a  point  lattice  composed  of  hypercubic 

unit  i_el  1  s . 

The  condition  expressed  in  equation  (3.3)  yields  important  insight  into 
the  judicious  choice  of  the  frequencies  «...  For  systems  with  large  dimensions, 

i  e.  a  large  number  n  of  rate  coefficients  (or  coupling  pa  ameters  in  general), 

1^1  approacties  m.  Thus,  the  cnoice  of  the  number  in  of  s  points  for  the  numerical 
computation  determines  the  value  of  u>  and  thus  guides  one  in  the  choice  of 
the  a.'s.  Since  one  would  expect  the  error  in  the  analysis  due  to  the  use  of 

ii  Leger  frequencies  u.  to  be  of  order  1/m,  i.e.  inversely  proportional  to  the  number 
of  jnit  cells,  it  is  evident  that  one  should  cnoose  a  large  value  o+  m  to  carry 

out  the  calculations.  Tnis  in  turn  implies  from  Eq.  (3.3)  a  large  [,j|. 

The  analysis  leading  to  Eq.  (3.3)  is  based  on  the  construction  of  unit 
cells  which  are  exact  hypercubes.  Since  |wj,  in  and  n  are  all  integers  it  may 
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not  be  possible  to  fulfill  condition  (3.3)  exactly  for  all  (arbitrary)  choices 
of  m  and  n.  It  is  easy  to  show  this .  Let  us,  for  instance,  square  both  bides 
of  Eq.  (3.3)  and  let  n  =  5.  We  are  then  required  to  find  an  integer  m  such  that 
we  can  express  a  sum  of  n  squares,  |w|^,  as  m^J.  Clearly,  this  is  possible, 
if  at  al lf  only  for  certain  special  values  of  in.  To  be  more  realistic,  we  should 
weaken  our  criteria  for  uniformity  of  the  point  lattice  by  stipulating  unit 
cells  which  are  as  close  as  possible  to  hypercubic  and  then  rewrite  Eq.  (3.3)  as 

Uli  (3.4) 

■  — ►  i  n 

M  =  n 

For  n  »  1,  tor  which  both  (3.3)  and  (3.4)  reduce  to  |u|  =  in,  it  should  be 
possible  to  obtain  a  more  nearly  exact  hypercubic  unit  cell. 

The  question  as  to  the  "optimum"  choice  of  the  w's  has  been  considered  by 
a  number  of  authors,  using  a  different  approach  from  tnat  presented  above,  in 
connection  with  the  general  problem  of  the  approximate  evaluation  of  multi¬ 
dimensional  integrals  via  discrete  summations^.  Korobov's  book  has  tables 
of  w^'s  for  a  given  number  of  points  m  (for  m  prime)  and  dimension  n  up  to 
n  =  10  These  taoles  are  reprinted  in  the  book  by  A.  H.  Stroud.  It  is  inter¬ 
esting  and  comforting  to  note  that  although  these  tables  were  computed  from 
completely  different  criteria  than  those  employed  by  us,  the  Korobov  2's  indeed 
generate  hypercubic  unit  cells  to  a  very  good  approximation  for  all  his  sets  of 
for  which  we  have  carried  out  the  appropriate  calculations. 

Korobov's  analysis  of  the  error  in  the  use  of  integer  frequencies  yields 
explicit  prescriptions  for  calculating  the  "optimum"  set  of  w.'s  as  given  in  his 
tables.  Our  approach  presented  above  does  not  yield  suen  an  explicit  algorithm. 
We  plan  to  develop  such  an  algorithm  and  then  compare  our  sets  of  ^  with  those 
of  Korobov  in  a  subsequent  publication. 
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frequency 

As  stated  above,  if  one  wishes  to  decrease  the  error  in  the  integer  A  anal ysis , 
one  should  use  a  large  number  m  of  s  points  and  thus  a  large  |u>|.  It  can 
readily  be  verified  from  Eq.  (2.6)  that  for  large  u.  one  needs  to  evaluate 
c(s)  for  a  larger  number  of  s  values  in  order  to  obtain  an  accurate  value  tor 
the  Fourier  coefficient  } .  This  requires  more  extensive  computer  calculations. 

A  reasonable  compromise  between  these  two  effects  needs  to  be  adopted. 

The  transformation  from  u  space  to  &  space  as  given  by  Lqs.  (2.3) 
will  also  effect  the  error  term  since  the  specific  transformation  which  is  chosen 
determines  the  rate  of  change  of  the  function  c(s)  as  a  function  of  s.  We  are 
faced  nere  with  an  interesting  problem  in  compensating  effects.  Either  the 
transformation  u.  -  u.(0)  is  singular  or  the  weight  function  p(u.;ai)  is  singu¬ 
lar  at  one  or  more  values  of  u..  For  the  chosen  transform  it  is  readily  veri¬ 
fied  from  Eqs .  (2.15)  and  (2.20)  that  u.  diverges  at  0  =  n/2  and  Q  =  3n/2,  but 
at  these  points  pfu^)  -  0.  This  same  effect  will  be  found  fop  aoy  transfor¬ 
mation  and  its  associated  weight  function.  Tnus  in  tne  regions  of  ft  space 
where  the  transformation  is  divergent  the  associated  weight  function  will 
always  compensate.  We  are  therefore  led  to  expect  that  the  choice  of  the 
transformation  function  will  not  significantly  effect  the  final  numerical  results. 
This  is  born  out  by  the  data  presented  in  paper  II. 
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